Build Classifier

Check a priori assignments and assign genetic clusters For the first step, will run successive k-means clustering after data transformation using PCA

Nice, k-means on the data finds two clusters (based on BIC), now let’s check to see if this fits with the assigned phenotypes of the adult fish.

# it appears the best cluster number is 2 so move on
grp <- find.clusters(adults, max.n.clust=10, n.pca=310, n.clust = 2)
table(pop(adults), grp$grp)

##         
##           1  2
##   Summer  0 41
##   Winter 40  0

K-means perfectly assigns genetic clusters to phenotypes!

Build the DAPC Classifier
Now that we did a little quality control to confirm that clusters in the PCA transformed data conform to the assigned phenotypes, we will conduct our DAPC.

First we run the dapc with all pcs retained to choose the number of pcs we need to retain info.

Seems like a small number of pcs is good, but let check how many are needed specifically for group assignment using the a-score (succesful reassignment probability after permuting group memberships)

#but just to check, lets run the a-score optimizer of adegenet
invisible(dapc_highpc <- dapc(adults, adults$pop, n.cluster=2, n.pca = 50))
optim.a.score(dapc_highpc, smart=FALSE)

This suggests that a single pc will adequately discriminate the groups without overfitting. Since we’re not concerned about overfitting (applying the trained classifier to new data), we’ll keep a few more than one PCs

#optim a score is just 1 pc
dapc_adult <- dapc(adults, adults$pop, n.pca = 5, n.da = 1)
scatter.dapc(dapc_adult, legend=TRUE)

Here we can see the DAPC separates known run timings on the basis of a single discriminant axis based on 5PCs

Examine the loadings
Yet another sanity check: Now we can check the loadings of SNPs on the PCs and DAs. Should return the known run history alleles. In addition to a sanity check this might be interesting to see if some chr28 SNPs are more informative than others for this population (Rogue)

marker_loadings <- loadingplot(dapc_adult$var.contr, axis=1,thres=.01, lab.jitter=1)

names(marker_loadings$var.values)

##  [1] "Chr28_11625241.1"  "Chr28_11625241.3"  "Chr28_11658853.1" 
##  [4] "Chr28_11658853.2"  "Chr28_11667578.4"  "Chr28_11667578.2" 
##  [7] "Chr28_11671116.2"  "Chr28_11671116.4"  "Chr28_11676622.4" 
## [10] "Chr28_11676622.3"  "Chr28_11683204.3"  "Chr28_11683204.4" 
## [13] "Chr28_11702210.4"  "Chr28_11702210.3"  "Omy_GREB1_05.4"   
## [16] "Omy_GREB1_05.3"    "Omy_RAD15709-53.1" "Omy_RAD15709-53.3"
## [19] "Omy_RAD47080-54.1" "Omy_RAD47080-54.3" "Omy_RAD52458-17.2"
## [22] "Omy_RAD52458-17.1"

Now, lets check if these are labeled as adaptive/run timing SNPs in the metadata file for the GTseq markers

#opened the "Omy GTseq panel updated Jun-2020.xlsx" file and converted into format for reading into R (pulled just the "Assay" and "Presumed type")
markersummary <- read_tsv("genotype_data/gtseq_marker_summary.txt")

# adegenet saves the genotype value as part of the SNP name, after a ".", strip this off first
clean_names <- str_remove_all(names(marker_loadings$var.values), ".[0-9]$")

markersummary[(markersummary$Assay %in% clean_names),]

missing info some of the marker names in the gtseq data are not in the marker summary spreadsheet I received, find the correct version

So far, only 4 of the 11 sites with high loading have matches in the gtseq summary file, but yes, all are annotated as run-time related

Classification

Now that we have a classifier, we’ll apply to to the half pounders.

pred.half<- predict.dapc(dapc_adult, newdata=half_pound)

# now lets get the individual data into a single df and plot 
ind.scores_adult <- as.data.frame(cbind(dapc_adult$ind.coord, as.character(dapc_adult$grp)))
ind.scores_half <- as.data.frame(cbind(pred.half$ind.scores, rep("halfpounder", times=length(pred.half$ind.scores))))
ind.scores <- bind_rows(ind.scores_adult, ind.scores_half)
ind.scores$LD1 <- as.numeric(ind.scores$LD1) #this is why you always use tidyR for datawrangling...

#now lets 

#plot
ggplot(data=ind.scores)+geom_density(aes(LD1, fill=V2, color=V2), alpha = 0.2)+scale_color_viridis_d()+scale_fill_viridis_d()+theme_classic()

ggplot(data=ind.scores)+geom_histogram(aes(LD1, fill=V2, color=V2), alpha = 0.2)+scale_color_viridis_d()+scale_fill_viridis_d()+theme_classic()

A single PC
Since assignment success of randomized groups was actually maximized ith a single pc summarizing differences among run types (see a-score analysis above). Let repeat the analysis with a single PC to see if we can improve assignment.

dapc_1pc <- dapc(adults, n.pca = 1, n.da=1)

pred.half.1pc<- predict.dapc(dapc_1pc, newdata=half_pound)

# now lets get the individual data into a single df and plot 
ind.scores_adult1 <- as.data.frame(cbind(dapc_1pc$ind.coord, as.character(dapc_1pc$grp)))
ind.scores_half1 <- as.data.frame(cbind(pred.half.1pc$ind.scores, rep("halfpounder", times=length(pred.half.1pc$ind.scores))))
ind.scores1 <- bind_rows(ind.scores_adult1, ind.scores_half1)
ind.scores1$LD1 <- as.numeric(ind.scores1$LD1) #this is why you always use tidyR for datawrangling...

#now lets 

#plot
ggplot(data=ind.scores1)+geom_density(aes(LD1, fill=V2, color=V2), alpha = 0.2)+scale_color_viridis_d()+scale_fill_viridis_d()+theme_classic()

ggplot(data=ind.scores1)+geom_histogram(aes(LD1, fill=V2, color=V2), alpha = 0.2)+scale_color_viridis_d()+scale_fill_viridis_d()+theme_classic()

About the same.

Rough Results Summary

Half pounders can not be reliably discriminated into winter or summer run timing on the basis of SNPs in the GTseq panel. They (loosely) fall into three groups: winter, summer and intermediate.

Formal Assignment
Given the distribution of known run timing individuals along the discriminant axis, we’ll (arbitrarily) define assignment of half pounders into run types as fish that fall into 95% credible intervals of known winter/summer fish

CIs <- ind.scores %>%
  group_by(V2) %>%
  summarise(loCI = quantile(LD1, probs = 0.025),
            hiCI = quantile(LD1, probs = 0.975))

## `summarise()` ungrouping output (override with `.groups` argument)

#number of half pounders that fall in the 95% credible interval for winter fish assignment
ind.scores %>%
  filter(V2 == "halfpounder") %>%
  summarise(winter_assigned = sum((LD1 < CIs$hiCI[3] & LD1 > CIs$loCI[3])), summer_assigned = sum((LD1 < CIs$hiCI[2] & LD1 > CIs$loCI[2])), unassigned = sum((LD1 < CIs$loCI[3] & LD1 > CIs$hiCI[2])) )

We should also see if using a single PC improves assignment:

CIs1 <- ind.scores1 %>%
  group_by(V2) %>%
  summarise(loCI = quantile(LD1, probs = 0.025),
            hiCI = quantile(LD1, probs = 0.975))

## `summarise()` ungrouping output (override with `.groups` argument)

#number of half pounders that fall in the 95% credible interval for winter fish assignment
ind.scores1 %>%
  filter(V2 == "halfpounder") %>%
  summarise(winter_assigned = sum((LD1 < CIs1$hiCI[3] & LD1 > CIs1$loCI[3])), summer_assigned = sum((LD1 < CIs1$hiCI[2] & LD1 > CIs1$loCI[2])), unassigned = sum((LD1 < CIs1$loCI[3] & LD1 > CIs1$hiCI[2])) )

Yes a single PC is more inclusive when assigning individuals to groups